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ABSTRACT 

John H. Bearden, Master of Science, 1977 
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ABSTRACT 


High Reynolds number, Incompressible flow has posed many 
problems, some seemingly insurmountable, to researchers investi- 
gating this phenomenon. It will be the purpose of this investigation 
to consider what has been done in this field, and use this as a 
basis to attempt to overcome these problems and achieve a simulation 
of incompressible high Reynolds number flow. 

A choice was made to use the stream function-vorticlty form 
of the Navier-Stokes equations. The coordinate system used was 
a body fitted coordinate system with a "U" shaped outer boundary. 

A modified version of a numerical solution computer program was 
used to actually solve the Navier-Stokes equations around the body 
used in the research. 

The particular body used in the Investigation was an NACA- 
0C18 airfoil. There was only moderate success, though, in this 
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simulation of flow at a Reynolds number of 1,000,000 and body angle- 


of-attack of zero. 
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CHAPTER I 
INTRODUCTION 

Since modern computers made the numerical solution of the 
Navier-Stokes equations practical , the goal of researchers in the 
Incompressible flow regime has been to simulate flow at high Reynolds 
numbers. The problems facing the researchers have proved to be 
excessive computer storage requirements as well as tremendous amounts 
of computer time needed. The computer storage requiremei^t stems 
from the need for having very closely spaced coordinate lines in 
the field over which the solution is made because of the high Reynolds 
number flow. The necessity of using great quantities of computer 
time is caused by the use of extremely small time steps in the solu- 
tion to maintain stability. Presented here is an approach to over- 
come at least part of these problems, and in doing so, achieve an 
accurate simulation of high Reynolds number flow about an airfoil. 

There have already been computer programs written that solve 
the Navier-Stokes equations about arbitrary bodies that work well at 
lov7 Reynolds numbers, so one of titese was used as a basis for this 
work. Also, a "U" shaped coordinate system has been developed for 
use at high Reynolds numbers which was used. The difficulties faced 
in the research were the conversion of the computer program that solved 
the Navier-Stokes equations to a form compatible with the "U" shaped 
coordinate system, the experimental evaluation of the various accel- 
eration parameters, and the determination of the body vorticity such 
that the body velocity would be zero. The work herein was done in 
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the stream function-vorticity form of the Navier-Stokes equations. 

A very thorough discussion of the actual computational aspects 

and the computer code is contained in Reference [ l] ; this work by 

Thames has proved to be the backbone of this research. Also of 

possible interest to the reader is work by R. N. Reddy. This work 

by Reddy is also in the high Reynolds number range, but was done in 

[ 2 ] 

the integro-dif ferential form. Work has also been done on sub- 
merged hydrofoils by S. P. Shanks. Shanks’s research produced a 

computer code to simulate the flow around a submerged hydrofoil, 

[31 

both with and without bottom effects. 
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CHAPTER II 

/mLYTIC DEVELOPMENT 

Due to the nature of this research, it might be advantageous to 
take a brief look at the derivation of the stream function-vorticity 
form of the Navier- Stokes equations and the transformation of these 
equations into the form for use in the transformed plane. A 
limited discussion of the related boundary conditions might also 
prove helpful. 

In cartesian coordinates, the Naviar-Stokes equations for non- 
steady, incompressible flow are: 

u + V = 0, [ 1-la ] 

X y 

U(. + u Ujj + V Uy = ~-^x [ 1,1b ] 

v+uv +VV = -ip + vV^v, [ l*lc ] 

t X y p y 

where u and v are the velocity components in the x and y directions, 
respectively, where the subscripts represent partial differentiation 
with respect to the variable used, and where p is the pressure. These 
equations are transformed to the stream function-vorticity form by 
first defining: 

u = tl(,v = -ii,a) = v - u; [1.2a)b,c] 

^y ^x’ X y’ 

then taking the curl of the Navier-Stokes equations, we have 

mt + i|'y ~ '^'x ^y ~ vV^co* f 1.3a ] 

= ~ Wy [ 1.3b ] 

With this Information, the next step will be to consider the 
transformation of the stream function-vorticity equations from the 
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form used in the physical plane to that used on the rectangular trans- 
formed plane. The equations are transformed using the following 
operators ; 


' ='n *1 ■ >'5 'n > ^ 


[2.4a] 


And upon utilizing these operators the following equations result: 

“t ^ ^ " '^5 “n ^ ^ ^ ^ 

( a 0)^^ - 2 p 0)^^ + Y + a (0^ + T (D^ ) / j2 tto, [ 2.5a ] 
a - 2 0 + Y + T 1//^ = - w, [ 2.5b ] 

where the functions m and i|) are doubly differentiable. The derivation 
of these expressions appears in Appendix A of Reference [ 1] . The 
coordinate system parameters a, 3, y» ^nd J are: 

“ = + yp ’ 

J = yr, - y^, 

3 = x^ Xr, + yg y^. 

T = r y| 

a 5 I y{ ( “ *55 - 2 B + Y !<„„ ) - 


( » yjl - 2 B yj, + Y y„J 1 / J 


T B ( < a 


2 3 y_ + Y y ^ “ 
Cn ^ nn 


y^ ( a Xg^ - 2 3 X^^ + Y x^^) ] / J 


Note that the Equations [2.5a,b] have been non-dlmensionalized 
witli respect to the free-stream velocity and the airfoil chord. The 
Reynold's number, R^, is also based on these quantities. 


5 


In the case under consideration, where the coordinate system 
remains fixed in time, the coordinate system parameters must only be 
calculated once since they remain constant throughout the solution. 
Also, if the coordinate system has no contraction of the coordinate 
lines about the body, the coordinate system parameters o and x 
disappear. 

The boundary conditions for the equations in the rectangular 
transformed plane are expressed in the following relations : 


'Kx,y,t) = i|;^(x,y,t) , x,y e F 2 

[2.7a] 

m(x,y,t) = oi^(x,y,t) , x,y e F^ 

[2.7b] 

'j'Cx,y,t) = iJjQ = constant, y,y e F^ 

[2.7c] 

• (x,y,t) - 0, [x,y] e Fi 

[2.7d] 


where Fj and V 2 are the body and the outer boundary line segments, 
respectively. Equation [2.7d] specifies conditions that, if met, 
guarantee the tangential velocity on the body surface to be zero. 

Since the normal component vanishes in the same manner, it can be 
assumed that the no-slip boundary conditions are met by satisfying 
this equation. 

For a much more indepth discussion on the material presented here 
the reader is to refer to Reference [1] . Here Thames covers the various 
aspects required in setting up the boundary conditions, as well as, the 
transformations used in manipulating the equations. 
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CHAPTER III 
THE COORDINATE SYSTEM 


The coordinate system used in this work is a body fitted 
coordinate system, that is, one that fits one coordinate line through 
every point specified on the body being considered. This type of 
coordinate line will be called an n-line. Another p-line is then 
fitted to the outer boundary and the remaining lines spaced throughout 
the field. A second type of lines is then used to subdivide the 
n-llnes, referred to here as £|-lines. Although the coordinate sys- 
tem created by this method is not necessarily orthogonal, there are 
no severe problems caused. On the contrary, the fact that the 
coordinate lines follow the body contour, thus eliminating any 
extrapolation, produces a profound advantage in the accuracy provided 
in the boundary layer. References [1], [4] , [5] , and [ 6 ] cover this 
topic in greater detail. 


Figure [1] shows such a coordinate system, in this case an 

airfoil in an elliptical outer boundary. Also, Figure [2] shov;s the 

rectangular transformed plane. This transformed plane is created by 

"opening up" the body fitted coordinate system from the cut. The 

line segment T in the body fitted coordinate system is analogous to 
* ^ 

in the transformed plane, the same correspondence also carries for 
* * * 
r 2 > as well as for and 1 ’^ and and Note that and 

are colinear, therefore each j^oint on 13 is equal to the corres- 
* 

ponding point on 

By starting with some initial guess, Figure [3], the body 
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fitted coordinate system is generated by solving two elliptic partial 
differential equations that have Dirichlet boundary conditions. 

These equations are: 


+ Cyy = P ( 5, n ) , 

[3.1a] 

+ Oyy = Q ( 5, n ) , 

[3.1b] 


with one coordinate being set constant on the body and outer boundary 
and the other varying monotonlcally around the body. Since the Initial 
guess specifies the values of x and y for each point, the 
dependent and independent variables in the equations must be inter- 
changed. Upon doing this, they take on the form of: 

a X -20 X +YX =-J^ x P(g, T))+x Q ( |,n ) [3. 2a] 

nn 5 n 

a yj-j- - 2 0 y + Y y = - J y- P ( €, h) + y Q ( ) [3.2b] 

55 5q pq 5 h 

This quasi-linear elliptic set of partial differential equations 
is much more complex than the linear equations [3.1], but the 
boundary conditions are specified at the body and outer boundary which 
are straight lines in the rectangular transformed plane. Also, the 
grid spacing in the transformed plane is unity. The inhomogeneoug 
terms in equations [3.2] , P(C,n) and Q (5>n)> are used in the 
contraction of the n~lines to the body, to a given line, or even to 
some specified point. These terms are derived from the sums of 
decaying exponentials. Further discussion of this topic is con- 
tained in Reference [&]• 

These equations are then solved over the rectangular transformed 
plane using central space differences and successive overrelaxatlon 
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iteration. 

With this as a basis, the coordinate system used in this research 
wi..! be considered. Instead of enclosing the body, in this case an 
airfoil, in an elliptical outer boundary, it is enclosed in a *'U" 
shaped outer boundary. Figure [4 ] . The reasons for doing this are 
the high contraction requirements and the necessity for good grid 
resolution in the wake. For high Reynolds number flow, the coordinate 
lines must be contracted very closely to the body to get enough lines 
in the boundary layer to man! tain stability in the solution. This 
creates problems, though, in that the lines that are contracted close 
to the body, in the system with an elliptical outer boundary, 
must bend around the sharp trailing edge. This bend brings the lines 
so close together that the separation becomes the order of magnitude 
of the roundoff error. But with the "U" shaped coordinate system, the 
n-lines "flow" off of the trailing edge completely eliminating this 
problem, Figure [6] . Also, the system with an elliptical outer 
boundary does not have sufficient lines in the wake area downstream 
of the body, and the lines spread out even more further downstream. 

The "U" shaped coordinate system also alleviates this problem 
because of the floxi? of the lines from the trailing edge. Note, 
also, that the lines have very 1/ttle spread dovmstream of the body. 
This gives excellent resolution in the wake. This type of system is 
created in the same manner as tl.e type vrith the elliptical outer 
boundary. The major difference is the way the rectangular trans- 
formed plane is set up, with the body and the reentrant segments 
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on the line and the infinity boundary occupying the other three 
sides. Referring to Figure [4] and Figure [5], some idea of the 
nature of the correspondence of the physical plane to the transformed 
plane may be gained. Again, in the physical plane is the same as 
the line segment in the transformed plane and p 2 corresponds to 
r*, to and so on. Since the coordinate lines F* and F* in the 
transformed plane represent the "cut" in the physical plane, they 
require special treatment when calculating the derivative across the 
"cut." Consider a point on that has coordinates (N,l), where 
1 < N LBYl, and where LBYl is the left hand coordinate of the body. 
It*s corresponding point on F* then has the coordinates [IMAX - (N-1) , 
1], where IMAX is the number of S“lines in the field. Using this, the 
derivative of ^ with respect to n becomes: 


'^'n ‘^1,2 ~ ’^IHAX-(N-1),2 

In the same manner the partial derivative of m with respect to 5 and n 
becomes : 


“Cn ~ \+l,2 ~ “iMAX- (N-2) , 2 ‘^IMAX-N,2 “ “n-1,2 
instead of Che form used in the field calculations which is: 


(i) ^ m.,, ~ to.,, , , + (0. , • 1 “ 1 -j.! 

1+1, j+1 1+1, j-1 1-1, J-1 i-1, j+1 

This holds true for all derivatives taken across the cut. 

The coordinate system was generated using the TOMCAT code in 


Reference [6]. The code permits tlie attraction of the lines to decay 
aft of the trailing edge thus permitting them to expand slightly, this 
feature was also used in resea ’-r h, The outer boundary was set up 5 
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chord lengths from the body except for the downstream dimension which 
was 10. The coordinate system has 10 coordinate lines contracted into 
the boundary layer, this attraction is explained in the appendix, 
Figure [6], 
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Perhaps the next step would be an explanation of the computer 
code used in this research to solve the Navier-Stokes equations over 
the rectangular transformed plane. This code was written by F. C. 
Thames, [1] and was chosen for the work because it uses the stream 
fiinction-vorticity form of the Navier-Stokes equations. 


The space derivatives of this solution are approximated using 
central differences and the vorticity time derivatives using a first- 
order backV'fard time difference. This implicit numerical method with 
SOR iteration is used to converge the space derivative variations in 
the flow solution. Then, with the field converged, the vorticity on 
the body is approximated with modified false position iteration. This 
method of iterating the body vorticity is designed to force the body 
tangential velocity to zero. ^ ^ ^ 


The numerical procedures used in accomplishing the solution and 
the finite difference form of the Navier-Stokes equations are now 
given attention. The finite dlflerenca grid, that Is the rectangular 
transformed plane, is pictured in Figure [7] for clarity. The grid 
pictured is one used in the *'U" shaped coordinate system. The dark- 
circles represent the body, on the J=1 line, and the outer boundary 
on the 1=1, I=IMAX, and J=JMAX lines. The points represented by the 

open circles are the re-entrant boundary segments and as pointed out 
earlier represent the same set of points in the physical plane. In 
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other words the point immediately down stream of the trailing edge 
in the physical plane is represented by the points I and II in the 
transformed plane, Figure [7 ] . Since the velocity at the outer boun- 
dary is the free stream velocity and the vorticlty is zero and since 
the body velocity is zero, the boundary conditions are specified at 
all points denoted by a darkened circle in Figure [7]. 


The procedure used in solving the difference equations over the 
rectangular transformed coordinate system is central differencing 
for the special derivatives and backwards differencing for the time 
derivative of the vorticlty. The difference equation for the vorti- 
city, after substituting the central difference and backwards dif- 
ference approximations, is: 

where the primed quantities represent difference expressions documented 

in Appendix D of Reference [1], Then solving for . yields; 

J 

* 

+ 


t 4.2 1 
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where; 

"=S)i,3 ^ 

® j/* 

(C5)j_j 5 44M'^^/4 
and 


The stream function equation is derived in the same manner and is: 

H- [<C3),,J - + 4ta - 2[(c;)^_j + (C5),_ ^/R„ 


The vorticity and stream function equations are then solved 
simultaneously over the computational grid by SOR iteration at each 
time step. Also it is necessary to calculate the vorticity iterate 
first in each sweep of the field if the two equations are coupled. 


As Thames points out, there are two observations that should be 
made about these difference equations. The first being the large 
quantities of computer storage needed for the implimentation of the 
solution. The reason being the nine field size arrays used to store 
the six coefficient arrays, the array, and the and ^ 


arrays . 
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Secondly, the fact that the velocity tangent to the body is explicitly 
enforced at zero alleviates instabilities in the solution. 

The changes made in this computer code were centered around the 
change from a coordinate system with the elliptical outer boundary to 
one that has the "U" shaped outer boundary. Recall from Chapter III 
that the re-entrant segments were moved from the 1=1 and I=IMAX lines 
in the transformed plane of the elliptical coordinate system to either 
side of the body on the J=1 line of the transformed plane of the "U" 
shaped coordinate system. Also the infinity boundary was moved to 
where it occupied the three remaining sides of the field of the "U" 
shaped coordinate system. Since no calculations are made on the 
infinity boundary, the original code swept the field from 1=1 to 
I=IMAX and then from J=2 to J=JMAX-1. That is every point is swept 
except the body and the infinity boundary at J=JMAX. But on the "U" 
shaped system, with the infinity boundary on the 1=1, I=IMAX, and the 
J=J>1AX lines, the field must be swept from 1=2 to I=IMAX-1 and J=2 
to J=JMAX-1, or in ohter words, all the interior points of the field. 
Then the points on the re-entrant boundary segments had to be calcu- 
lated. This was done by sweeping from 1=2 to the point immediately 
down stream of the body. Note that the J=-l points shown in Figure 
[8] correspond to the J=2 points of the other re-entrant segment. 

And, finally, the program was changed to sweep only the body points 
on the 1=1 lines and not the entire line. 


15 


The body vortlclty was calculated using Isreali's method for 
the first two changes in body vorticity. The algorithm is; 

But the method is slow and is therefore used only to make the first 
couple of changes in body vorticity. The subsequent iteration in body 
vorticity is done with the false position method; 


(k+1) (k) . 

■ “i.i - 


X. > X » 


t— - [- 


3n 


(n^) 


r — — 1 

Si /"-u 


(k) 


i,l 


8n 


(nj^) 


3n 


1,1 


1,1 


Because the boundary condition is Imposed on i|< and the vorticity, 

on 

(i), is calculated such that the body velocity is zero. R. L. Walker 
also used this method in an investigation of flow over a semi-infinite 
flat plate, Reference [8]. 
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CHAPTER V 
RESULTS 

The test case considered in this research was the flow about an 
NACA-0018 airfoil at a Reynolds number of 1,000,000, Detail of the 
coordinate system near the airfoil is shown in Figure [9], The entire 
field is shown in Figure [6], This case was run usini; a one-hundred 
step gradual start. The velocity was gradually increased by an in- 
crement of .01 from zero to free-stream velocity while holding the 
Reynolds number fixed at 1,000,000. The field vorticity acceleration 
parameter was 1.8, and the vorticity convergence tolerance, the stream 
function convergence tolerance, and the maximum body velocity were 
set at 0.0001, 0.00001, and 0.0001, respectively. 

There were some difficulties encountered in starting the solution. 
The primary difficulty was finding an accurate approximation for the 
trailing edge vorticity. The cause of the problem appeared to be a 
loss of symmetry of the vorticity between the top and the bottom of the 
airfoil, probably due to the extreme gradients encountered in the 
start. The number of Iterations required to converge the field was so 
large that the body velocity was converged to a value of .005 rather 
than .0001 during the start. Even with this change, it became necessary 
to ignore the value of velocity at the trailing edge points and the two 
points adjacent to the trailing edge. The vorticity at the trailing 
edge, then had to be set at zero rather than calculated because the 
velocity was not converged at the points used for extropolatlng the 
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vorcicity. This is valid for a symmetric flow. Finally, just before 
the solution reached free-stream velocity, the assymmetry caused the 
solution to diverge. To correct this, the points above and below the 
airfoil were averaged to make the field symmetrical. 

Once the solution reached free-stream velocity, the assytranecry 
was still present as shown in Figure [10] which pictures the stream- 
lines around the airfoil at time step 110. Since this assymmetrlc 
vorticity distribution caused the solution to diverge, it had to be 
somehow eliminated. Refer to Figure [11] and Figure [12] which show 
the velocity distribution around the airfoil. Note that the velocity 
leaving the upper surface is greater than that leaving the lower sur- 
face, even at time step 20. This could be caused by one of two things, 
first, the relative field vorticity iteration error was only converged 
to a value of .001, rather than to .0001, before the first change in 
body vorticity was made, thus leading to error before the solution is 
even started. Secondly, the field and body numerical iteration scheme 
sweeps from the trailing edge around the under side of the airfoil then 
over the top of the airfoil and back to the trailing edge. It is 
possible that this method of "sweeping" the field is perhaps biasing 
the velocity around the body, as a result of the large gradients at the 
start, leading to a higher velocity on the top of the airfoil. Figure 
[13] shows the pressure distribution on the airfoil at four times. Due 
to the above-mentioned numerical error, the assymmetric pressure dis- 
tribution becomes worse as time progresses. The pressure coefficient 
after the acceleration has stopped does, however, resemble the expected 
profile for this airfoil. The drag is about ice the experimental 
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value, but was still decreasing at the end of the run. The high drag 
and the linear tendency in the pressure distribution during the start 
are, of course, due to the acceleration. 

Two proposed changes to correct these problems are: (1) Converge 

the field vorticity iteration error more completely before attempting 
a change in body vorticity. (2) Sweep the field from the center to the 
outside edges, i.e. , leading edge to trailing edge rather than from 
left to right, and start from a potential flow. 
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In this appendix, the procedure by which a specified number of 
coordinate lines can be automatically concentrated into a boundary layer 
of specified thickness is discussed. Consider the coordinate system 
generation equations (3.2) applied to the one-dimensional case of 
straight boundaries parallel to the x-axis. With n = constant on these 
boundaries, and the 5-lines being normal to the boundaries, we have 

= Y^^ = 0 and the x-equation is identically zero so that the 
coordinate equations reduce to 


vY + J2 qY == 0 

' nn n 


(A-l) 


or 


Y t2 

ja+ JI q . 0 

\ r 


(A. 2) 


This can be made a perfect differential by choosing the form of the 
control function Q to be 


Q(n) 


1 

7 ' f(n) 


where the minus sign has been introduced merely for convenience. 
(A. 2) becomes 


1 c t ^ 

V t 


(A. 3) 


Then 


(A. 4) 


*Personal communication from Dr. Thompson 
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which can be integrated to yield 

y(n) = + ^2 

The constants of integration may be evaluated from the boundary condi- 
tion: y(l) y(J) = yj so that 


y(h) = yj^ + (y^-y^^) 


f(n) - f(l) v 
4(J) - f(l)'' 


This equation should be solved for f(n) to yield 


(A. 6) 


f(n) - HI) ~ ^1 

f(J) - f(l) Yj - 


(A. 7) 


which, with arbitrary definition of f(l) and f(J) will yield the 
required f(n)j and hence the required Q(ii) via substitution in (A. 3), 
to produce a desired distribution y(n)- The evaluation of Q(n) may be 
done without actual evaluation of f(n), however, by solving (A. 4) for 
y"/y’ and substituting into (A. 3) to pi-oduce 

Q(n) = - ^ (A. 8) 

j y 

Now a number of smooth functions for y(n), such as exponentials, 
logarithmic functions, hyperbolic function, etc,, may be found which 
will concentrate lines near y^ with a spread out to y^. However, since 
the boundary layer thickness at high Reynolds number is only a very 
small fraction of the distance to outer boundary of the computational 
field, such smooth functions cannot allow the lines to spread rapidly 
enough outside of the boundary layer. The result is that nearly all 
of the lines in the field will be within a few boundary layer thick- 
nesses of the body, with a great gap near the outer boundary. 
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Therefore, a composite function was used for y(n)i formed by 
joining a logarithmic function to a quartic polynomial near the edge of 
the boundary layer. This function was constructed as follows: assume 

that it is “desired to space the lines in the boundary layer such that 
the change in velocity from each at the next is the same. Let the 
velocity profile in the boundary layer be approximated by the expo- 
nential 

u(y) = 1 - e“®^ (A. 9) 

Let the edge of the boundary layer be defined by 

u = 0.99 at y = 6 

Then the decay factor c will be given by 

c = - •! In (0.01) (A. 10) 

Now solve (A. 9) for y(u): 

y(u) = - In (1-u) (A. 11) 

In order to achieve the same velocity change from each line to the 
next, take u = 0.99 ( ^ " v ) where rit. is the line at the edge of the 
boundary layer. Substitution is (A. 11) then yields 

y(n) = - - In [1 - 0.99 (-^) ] 1 < n < n. (A. 12) 

c *^5 ° 

Let this logarithmic function be joined to a quartic polynomial 
at some line inside or at the outer edge of the boundary layer. Thus 
with the function at n=N, the polynomial is of the form 

y(n) =y'(N) [n-N] +|y"(N) [n~N]2 

+ J y"(N) [n-N] 3 + a(n-N)‘‘ + y(N) 

N < n < J (A. 13) 
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Here y'(N) is functional notation, etc. The derivatives are determined 
by differentiation of (A. 12) with evaluation at n=N. The remaining 
coefficient "a" is used to satisfy the boundary condition at the outer 


boundary y(J) “ Yj* Thus 


a = 


Vj - y(N) - y'(N) [J-N] - I y"(N) [J-N]^ - i y"(N) [J-N]3 


(J-N)' 


(A.14) 


Note that the junction to the polynomial need not occur at the edge of 
the boundary layer, but anywhere inside it. It was found advantageous 
to place the junction two or three lines inside the boundary layer. 

Thus if the boundary layer thickness, 6, and the number of lines 
therein, , are specified, along with the distance to the outer 
boundary, y^, and the total number of lines J, and the function line N, 
the control function Q(n) can be evaluated from 

0.99 

Y ^(5 

Q(n) = - n = 1,2, — , N £ (a, 15) 

1 - 0,99 ,-) 

n -1 


Q(n) = - 


y"(H) + y"(N) [n-N] + 12a[il-N]2 


y' (N) + y”(N) fn-N] + ■- y”(N) [n-N]^ + 4a[n-N]^ 


n = N, N + 1, — , J 


(A. 16) 


with the required derivatives given by 


y^^^ (N) 


(k-1)! .O^^.k 
c '■ W-1 ^ 


[1 - 0.99 


k = 1,2,3 


(A. 17) 
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and y(N) by 

y(N) = - ln[l - 0.99 (A.18) 

c Hg-l 

Although this analysis is developed for the one-dimensional case of a 

flat boundary, much the same results will be achieved by its use with 

curved boundaries since curvature tends to affect both the boundary 

layer thickness and the line control in the same way. Thus convex 

curvature thins the boundary layer but also causes the lines to 

concentrate to a greater degree near the boundary. 

In the present work, the boundary layer thickness was taken as 
5 

(S = 7~R where R is the chord Reynolds number, and 10 lines were placed 
therein (n^ = 10) with the junction to the polynomial at line 7 (N = 7). 
There were 31 lines in the field, and the outer boundary was at 5 
(J = 31, yj = 5). 

Additional coordinate system control, in the form of ?-line 
attraction was used to pool the C“lines in the wake nearer the trailing 
edge. This attraction was of the exponential type used in Thames [1] 
and in the original TOMCAT code [6], except that in order to be 
compatible with the boundary layer attraction function, the attraction 
was calculated for 

Q(n) = - (, S(n) (A. 19) 

J 

where S(n) corresponds to the Q(ii) of [1] and [6]. The attraction was 
point attraction to the trailing edge, with amplitude of 0.7 on the 
bottom of the trailing edge and -0.7 on the top. The decay factor was 
0.1 and the feature of attraction to the convex side and repulsion to 
the concave provided for in the TOMCAT code was activated. 
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Figure 1. Coordinate System - Elliptical Outer Boundary 






Figure 3. Initial Guess for Coordinate Generation 
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Figure 4, Coordinate System - "U" Sha*;ed Outer Boundary 
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Figure 8 . Section of the Computational Grid at the Cut 





Figure 10. Streamlines (t ■ 1.1) 
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Figure 13- Pressure 
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